Disentangling the critical signatures of neural activity

The critical brain hypothesis has emerged as an attractive framework to understand neuronal activity, but it is still widely debated. In this work, we analyze data from a multi-electrodes array in the rat’s cortex and we find that power-law neuronal avalanches satisfying the crackling-noise relation coexist with spatial correlations that display typical features of critical systems. In order to shed a light on the underlying mechanisms at the origin of these signatures of criticality, we introduce a paradigmatic framework with a common stochastic modulation and pairwise linear interactions inferred from our data. We show that in such models power-law avalanches that satisfy the crackling-noise relation emerge as a consequence of the extrinsic modulation, whereas scale-free correlations are solely determined by internal interactions. Moreover, this disentangling is fully captured by the mutual information in the system. Finally, we show that analogous power-law avalanches are found in more realistic models of neural activity as well, suggesting that extrinsic modulation might be a broad mechanism for their generation.

T 0 s(t, T )dt. At the critical point we expect that plots of t/T versus s(t, T )T 1−δ for different T will collapse onto the same universal scaling function [1].
Thus, finding the exponent for which the goodness of the collapse is higher provides another way to estimate δ. For testing the avalanche shape collapse, we used the methodology introduced in [2]. To determine the quality of the collapse, the averaged and rescaled avalanche profiles of different lifetimes F (t/T ) = T 1−δ s(t/T, T ) are first linearly interpolated at 1000 points along the scaled duration. The variance across the different F (t/T ) is calculated at each interpolated point, and the shape collapse error ϵ(δ) is then defined as the mean variance divided by the squared span of the avalanche shapes, where the span equals the maximum minus the minimum value of all rescaled avalanche profiles. In the presented analysis, avalanche shapes of T > 10 bins with at least 10 samples were used.
The collapse has been tested on the extrinsic model in the case of low D * , and the results are plotted in Figure S1. We find that the exponent that minimizes ϵ(γ) turns out to be ≈ 1.33, close to the estimates of δ found through the linear fit of average size given duration and through the prediction of the crackling-noise relation in the main text. Again, it is also close to the apparently super-universal exponent found in [3] and in [4].
FIG. S1. Collapse of the average profile of avalanches of varying duration in the extrinsic model, for the low D * regime. (a) Profile of the avalanches before the rescaling. (b) If we rescale with an exponent δ ≈ 1.33, which is remarkably close to the one found in the main text through the crackling-noise relation, we obtain an optimal collapse onto the same scaling function.

TESTING THE PREDICTIONS OF THE EXTRINSIC MODEL
Here, we show that simulations of the extrinsic model described by an Ornstein-Uhlenbeck process agree with the analytical results presented in the main text. Whenever not specified, we assume that the parameters of the model are given by D * = 0.3, θ = 1, γ D = 10 together with γ i = 0.1, γ j = 0.5. Thus, we are in the limit of timescale separation considered in the main text.    S4. Comparison between the analytical expressions of the joint probability distribution p(v i , v j ) given by Eq. (S1) and the one obtained from with 10 3 simulations. We also show the corresponding results for the factorized probabilities. The blue line corresponds to the analytical expression of p(v i , v j ). The corresponding dots represent the histogram of the distribution obtained from 10 3 simulations, and the semitransparent filled areas represent one standard deviation from this estimate. Similarly, the gray dashed line represent the analytical expression of p(v i )p(v j ). (a) Section along the v j direction for small v j , so that we are looking at the tails of the distribution. Even though the estimate along the tails is noisy, we clearly see that the estimate from the simulations lies along the analytical prediction. (b) As before, but for higher v j . (c-d) As before, but with values of v j close to zero so we look at the bulk of the distribution near its peak. Even though joint probability and its factorization now are more similar, once more the estimate from the simulation match the analytical expression p(v i , v j ).

Let us recall that
is the stationary joint probability of the model. In Figure S2a we check that the stationary distribution of D is indeed the one presented in the main text and in Figure S2b that the stationary distribution of a single v i corresponds to the analytical expression we derived. If we compare this distribution to a standard distribution of an Ornstein-Uhlenbeck process with a diffusion coefficient equal to the mean ⟨D(t)⟩ we immediately see that the distribution of our model is considerably more peaked around zero and displays longer tails. Indeed, one expects that due to the fact that D * < ⟨D(t)⟩ the system tends to wander more easily close to zero, especially in the time windows where the diffusion coefficient is constant and equal to D * . At the same time, the fact that D(t) can change in time favors the presence of values of v that are larger in absolute value, which is the mechanism at the origin of the bursty behavior seen in the main text.
In Figure S3 we look instead at the properties of the joint probability distribution p(v i , v j ). The most natural quantity to compare this distribution with is its factorization p(v i )p(v j ), which is equivalent to ignoring the feedback effects between v i and v j due to the shared extrinsic modulation of D(t). Since we are setting D * = 0.5, we expect that these effects are going to be particularly relevant for the dynamics of the model. In particular, in Figure S3c-d we see that the most important differences between the two occur in the tails of the two-dimensional distribution, with the joint distribution typically showing dramatically longer tails. This translates to the fact that far-from-zero values of the two variables can occur more easily at the same time. The situation is completely reversed when we increase D * . In Figure S3e-f we see that for D * = 8 the joint probability distribution and the factorized distribution are almost indistinguishable. Hence, this example shows explicitly that if D * is high enough the dependence induced by the extrinsic modulation vanishes.
Let us keep focusing on the case D * = 0.5 for the time being. In Figure S4a-d we compare the one-dimensional sections of the analytical expression of the joint probability distribution with the results of 10 3 simulations of the model, together with the sections of the factorized distribution. The joint distribution estimated from the simulation matches particularly well the analytical prediction. Once more, and perhaps more clearly, in Figure S4a-b we see the stark difference that emerges along the tails between the joint probability distribution and its factorization. Interestingly, panel (c) and panel (d) show that the situation in the bulk of the distribution is reversed with respect to the tails and now the joint probability distribution is more peaked with respect to is factorization, albeit only slightly. That is, the modulation in the low D * regime favors both large values of v i and v j and values very close to zero.
Overall, this brief analysis showed us how the bursty behavior that we see for small values of D * emerges from the underlying probability distributions, which in turn emerge from the simple marginalization that occurs in Eq. (S1). Similar arguments, albeit impractical, could be carried out for the probability distributions beyond the two-point ones. In a sense, one could argue that the fundamental properties of the model are inherited from the fact that there are some unobserved physical quantities, and these are the quantities that drive the global response of the single variables.
For the specific case of a double Ornstein-Uhlenbeck process we chose, the net effect of the marginalization is the widening of the tails of both the one-point p(v) and the two-point p(v i , v j ) probability distributions when D * is small enough. As we increase D * , this effect becomes less and less important until it is completely negligible. In this sense, we can effectively think of D * as a control parameter that changes the qualitative behavior of the system. Most importantly, the fact that the tails of the joint probability distribution are wider when D * is small reflects dynamically in the emergence of a non trivial coordination between the variables, from which in turn power-law avalanches emerge.

CORRELATIONS IN THE INTERACTING MODEL
In Figure S5 we compare the correlation matrix obtained from the data, the interaction matrixÃ ij obtained by solving the inverse problem described in the main text, and the correlation matrix obtained from simulations of the resulting model.  Figure S6 shows avalanches' distributions in the extrinsic model as D * is increased. It is possible to see that, while the exponents converge to a fixed value for D * sufficiently low, as D * is increased they gradually increase in absolute value, and the distribution tends to an exponential when D * > 5. That is, the shift from exponential and power-law avalanches is smooth. Figure S7 shows instead the finite-size scaling behavior of avalanches in the extrinsic model. The avalanches obtained do not present a clear cutoff. Yet, it is evident that when increasing the number of units the values of the maximum avalanches size and duration increase, while the exponent of the distribution remains the same, displaying a property typical of true power laws. This is perhaps unsurprising, since the units are conditionally independent to begin with. In [5] (SI) we have studied the avalanche finite-size scaling in LFPs from the same set of experiments, and they indeed reproduce the expected finite-size scaling features, as already found in the literature [6]. Avalanches are now fitted with an exponential distribution. (h) The average avalanche size as a function of the duration scales with an exponent that, as σ (h) decreases, becomes closer to the trivial one δ fit ≈ 1.

A CRITICAL STOCHASTIC WILSON-COWAN MODEL AS EXTERNAL MODULATION
In Figure S8 we analyze some Wilson-Cowan units that receive the same input by another Wilson-Cowan population that is in a critical state. Indeed, in this case we set α (h) = ω = 0.1, and we are exactly at the critical point of the stochastic Wilson-Cowan model as shown in [7]. Here, avalanches in principle will be present for all noise amplitudes (i.e., system sizes). However note that in Figure S8 we set h (h) = 10 −3 , so that the dynamics slightly deviates from the critical point defined for h = 0. For this reason, avalanches disappear if we reduce the noise to a sufficiently low value (see Figure S8e-h). Notably, the noise amplitude has to be dramatically reduced in order not to see avalanches, at σ (h) = 5 × 10 −5 , corresponding to very large system sizes (K = L ≈ 4 × 10 7 , with K and L respectively the number of excitatory and inhibitory neurons [8] corresponding to the neural population).

AVALANCHES STATISTICS ACROSS RATS
In Figure S9 we report the avalanches statistics from the other rats that we analyzed. The avalanches statistics is computed by considering all available the 20 trials of basal activity for each rat, that are 7.22s long. Inter-rat variability is present with respect to avalanche exponents, and is expected as found in previous experiments [3,9,10]. Moreover, a theoretical explanation for the difference in avalanche exponents has been also recently proposed as signature of quasi-criticality [11]. However, the fundamental point is that the crackling-noise relation is always verified compatibly with the experimental errors, a feature that is usually considered a hallmark of criticality [3].  (1a-4a) The distribution of the avalanches' sizes is consistently a power-law, with an exponent that slightly depends on the single rat. (1b-4b) The avalanche durations are once more power-law distributed in all rats with some variability in the exponents, even though the range accessible with the experimental setup only covers two decades. (1c-4c) The crackling-noise relation, however, is consistently satisfied in each rat.